"""
    GSIIntegration

Integration module for DRP-4DVar with GSI atmospheric data assimilation components.

This module provides interfaces that allow the DRP-4DVar algorithm to work
with realistic atmospheric observation operators, background error covariance
models, and atmospheric forecast models used in operational weather prediction.

# Key Components

- `GSIObservationInterface`: Adapters for GSI observation operators
- `GSIBackgroundErrorInterface`: Integration with GSI background error covariance
- `GSIModelInterface`: Model operator integration for 4D-Var time evolution
- `AtmosphericTestCase`: Realistic test case generators
- `OperationalDRP4DVar`: Operational-ready DRP-4DVar configuration

# Features

- Realistic atmospheric observation types (radiosondes, aircraft, satellite)
- Operational background error covariance modeling
- WRF/GFS model integration
- Regional and global domain support
- Multi-scale localization
- Adaptive ensemble inflation

"""
module GSIIntegration

using LinearAlgebra
using SparseArrays
using Printf
using Statistics
using Random

# Import from parent modules (conditionally to avoid circular dependencies)
using ..FourDVar

# Note: These modules would be imported in a real GSI system
# For now, we'll use simplified implementations or skip dependencies
# using ..GridOperations

# Background-error handling relies on the core GSICoreAnalysis modules.  Use
# triple-dot relative imports because GSIIntegration lives under FourDVar.
using ...BackgroundError
using ...ControlVectors
using ...GSICoreAnalysis: AnalysisConfig, AbstractAnalysisConfig

export GSIObservationInterface, GSIBackgroundErrorInterface, GSIModelInterface
export AtmosphericTestCase, OperationalDRP4DVar
export create_realistic_observations, create_gsi_background_error
export apply_background_covariance, apply_background_covariance_inverse
export sample_background_covariance
export create_atmospheric_test_case, run_operational_drp4dvar
export add_temperature_observations!, add_wind_observations!, add_pressure_observations!, add_model_operators!
export create_observation_operators

# =============================================================================
# GSI Observation Interface
# =============================================================================

"""
    GSIObservationInterface

Interface for integrating GSI observation operators with DRP-4DVar.

This structure wraps GSI observation operators to provide the function
signature expected by the DRP-4DVar algorithm while maintaining
compatibility with realistic atmospheric observations.
"""
struct GSIObservationInterface
    # Observation operators by type and time
    temperature_ops::Dict{Int, Any}
    humidity_ops::Dict{Int, Any}
    wind_ops::Dict{Int, Any}
    pressure_ops::Dict{Int, Any}
    satellite_ops::Dict{Int, Any}
    
    # Grid information
    grid::Any
    
    # Observation metadata
    obs_locations::Matrix{Float64}
    obs_types::Vector{Int}
    obs_times::Vector{Int}
    
    function GSIObservationInterface(grid; 
                                    obs_locations=Matrix{Float64}(undef, 0, 3),
                                    obs_types=Int[],
                                    obs_times=Int[])
        new(
            Dict{Int, Any}(), Dict{Int, Any}(), Dict{Int, Any}(),
            Dict{Int, Any}(), Dict{Int, Any}(),
            grid, obs_locations, obs_types, obs_times
        )
    end
end

"""
    add_temperature_observations!(interface, locations, times, operator_config)

Add temperature observations to the GSI observation interface.
"""
function add_temperature_observations!(interface::GSIObservationInterface,
                                     locations::Matrix{Float64},
                                     times::Vector{Int},
                                     operator_config=nothing)
    
    nx, ny, nz = _grid_dimensions(interface.grid)
    
    for time_idx in times
        interface.temperature_ops[time_idx] = function(state_vector)
            n_3d = nx * ny * nz
            t_field = reshape(state_vector[(2*n_3d+1):(3*n_3d)], nx, ny, nz)
            return _sample_field(t_field, locations, nx, ny, nz)
        end
    end
    
    # Update metadata
    n_new_obs = size(locations, 1) * length(times)
    current_n = size(interface.obs_locations, 1)
    new_locations = vcat(interface.obs_locations, repeat(locations, length(times), 1))
    new_types = vcat(interface.obs_types, fill(220, n_new_obs))  # Temperature obs type
    new_times = vcat(interface.obs_times, repeat(times, inner=size(locations, 1)))
    
    # Update interface (this would need mutable struct in practice)
    @warn "Note: GSIObservationInterface should be mutable for this operation"
    
    return nothing
end

"""
    add_wind_observations!(interface, locations, times, component, operator_config)

Add wind component observations to the GSI observation interface.
"""
function add_wind_observations!(interface::GSIObservationInterface,
                               locations::Matrix{Float64},
                               times::Vector{Int},
                               component::Symbol=:u,
                               operator_config=nothing)
    
    field_offset = component == :u ? 0 : 1  # u is first field, v is second
    nx, ny, nz = _grid_dimensions(interface.grid)
    
    for time_idx in times
        interface.wind_ops[time_idx] = function(state_vector)
            n_3d = nx * ny * nz
            slice_start = field_offset * n_3d + 1
            slice_end = (field_offset + 1) * n_3d
            wind_field = reshape(state_vector[slice_start:slice_end], nx, ny, nz)
            return _sample_field(wind_field, locations, nx, ny, nz)
        end
    end
    
    return nothing
end

"""
    add_pressure_observations!(interface, locations, times, operator_config)

Add surface pressure observations to the GSI observation interface.
"""
function add_pressure_observations!(interface::GSIObservationInterface,
                                   locations::Matrix{Float64},
                                   times::Vector{Int},
                                   operator_config=nothing)

    nx, ny, nz = _grid_dimensions(interface.grid)
    n_3d = nx * ny * nz
    n_2d = nx * ny

    for time_idx in times
        interface.pressure_ops[time_idx] = function(state_vector)
            ps_start = 4 * n_3d + 1
            ps_end = ps_start + n_2d - 1
            ps_field = reshape(state_vector[ps_start:ps_end], nx, ny)
            return _sample_field(ps_field, locations, nx, ny, 1)
        end
    end

    return nothing
end

function _grid_dimensions(grid)
    nx = hasproperty(grid, :nx) ? getproperty(grid, :nx) : size(grid, 1)
    ny = hasproperty(grid, :ny) ? getproperty(grid, :ny) : size(grid, 2)
    nz = hasproperty(grid, :nsig) ? getproperty(grid, :nsig) : (hasproperty(grid, :nz) ? getproperty(grid, :nz) : 1)
    return nx, ny, nz
end

function _sample_field(field::AbstractArray{T,3}, locations::Matrix{Float64}, nx::Int, ny::Int, nz::Int) where T
    n_obs = size(locations, 1)
    samples = Vector{T}(undef, n_obs)
    for obs_idx in 1:n_obs
        i = clamp(Int(round(locations[obs_idx, 1])), 1, nx)
        j = clamp(Int(round(locations[obs_idx, 2])), 1, ny)
        k = clamp(Int(round(locations[obs_idx, 3])), 1, nz)
        samples[obs_idx] = field[i, j, k]
    end
    return samples
end

function _sample_field(field::AbstractArray{T,2}, locations::Matrix{Float64}, nx::Int, ny::Int, _) where T
    n_obs = size(locations, 1)
    samples = Vector{T}(undef, n_obs)
    for obs_idx in 1:n_obs
        i = clamp(Int(round(locations[obs_idx, 1])), 1, nx)
        j = clamp(Int(round(locations[obs_idx, 2])), 1, ny)
        samples[obs_idx] = field[i, j]
    end
    return samples
end

"""
    create_observation_operators(interface) -> Dict{Int, Any}

Convert GSI observation interface to DRP-4DVar observation operators.
"""
function create_observation_operators(interface::GSIObservationInterface)
    operators = Dict{Int, Any}()
    
    # Combine all observation types for each time
    all_times = unique(vcat(
        collect(keys(interface.temperature_ops)),
        collect(keys(interface.humidity_ops)),
        collect(keys(interface.wind_ops)),
        collect(keys(interface.pressure_ops)),
        collect(keys(interface.satellite_ops))
    ))
    
    for time_idx in all_times
        # Create combined operator for this time
        operators[time_idx] = function(state_vector)
            observations = Float64[]
            
            # Apply each available operator
            for (op_dict, name) in [
                (interface.temperature_ops, "temperature"),
                (interface.humidity_ops, "humidity"), 
                (interface.wind_ops, "wind"),
                (interface.pressure_ops, "pressure"),
                (interface.satellite_ops, "satellite")
            ]
                if haskey(op_dict, time_idx)
                    obs_part = op_dict[time_idx](state_vector)
                    observations = vcat(observations, obs_part)
                end
            end
            
            return observations
        end
    end
    
    return operators
end

# =============================================================================
# GSI Background Error Interface
# =============================================================================

"""
    GSIBackgroundErrorInterface

Interface for integrating GSI background error covariance with DRP-4DVar.

This provides realistic background error covariance matrices based on
climatological statistics and ensemble-based error modeling used in
operational data assimilation systems.
"""
mutable struct GSIBackgroundErrorInterface
    # Background error configuration
    config::Any  # BackgroundError.BackgroundErrorConfig
    grid::Any
    
    # Error statistics
    error_variances::Dict{Symbol, AbstractArray{Float64}}
    correlation_lengths::Dict{Symbol, Float64}
    
    # Incremental-form background covariance operator
    analysis_config::AbstractAnalysisConfig
    background_operator::BackgroundError.AbstractBackgroundError
    
    function GSIBackgroundErrorInterface(config, grid, analysis_config,
                                         background_operator)
        variances = Dict{Symbol, AbstractArray{Float64}}()
        correlations = Dict{Symbol, Float64}()

        corr_config = _get_config_field(config, :correlation_lengths, nothing)
        correlations[:horizontal] = _get_config_value(corr_config, :horizontal, 100.0)
        correlations[:vertical] = _get_config_value(corr_config, :vertical, 1.5)

        new(config, grid, variances, correlations, analysis_config, background_operator)
    end
end

# -----------------------------------------------------------------------------
# Internal helpers for configuration extraction
# -----------------------------------------------------------------------------

function _get_config_field(config, key::Symbol, default)
    if config === nothing
        return default
    elseif config isa AbstractDict
        return get(config, key, default)
    elseif hasproperty(config, key)
        return getproperty(config, key)
    else
        return default
    end
end

function _get_config_value(config, key::Symbol, default)
    if config === nothing
        return default
    elseif config isa AbstractDict
        return get(config, key, default)
    elseif hasproperty(config, key)
        return getproperty(config, key)
    else
        return default
    end
end

"""
    apply_background_covariance(bg_interface, x)

Apply the implicit background error operator B to a state increment vector.
"""
function apply_background_covariance(bg_interface::GSIBackgroundErrorInterface,
                                     x::AbstractVector{<:Real})
    cv = ControlVector(bg_interface.analysis_config)
    @assert length(x) == length(cv.values) "State vector size mismatch with background operator"

    cv.values .= x
    result_cv = BackgroundError.apply_background_error(bg_interface.background_operator, cv)

    return copy(result_cv.values)
end

"""
    apply_background_covariance_inverse(bg_interface, x; regularization=1e-6)

Apply the inverse of the background error operator B⁻¹ to a state increment vector.
"""
function apply_background_covariance_inverse(bg_interface::GSIBackgroundErrorInterface,
                                             x::AbstractVector{<:Real};
                                             regularization::Real = 1e-6)
    cv = ControlVector(bg_interface.analysis_config)
    @assert length(x) == length(cv.values) "State vector size mismatch with background operator"

    cv.values .= x
    result_cv = BackgroundError.apply_background_error_inverse(
        bg_interface.background_operator, cv;
        regularization = float(regularization)
    )

    return copy(result_cv.values)
end

"""
    sample_background_covariance(bg_interface, ξ)

Project a white-noise control vector ξ through B^(1/2) to generate a correlated perturbation.
"""
function sample_background_covariance(bg_interface::GSIBackgroundErrorInterface,
                                      ξ::AbstractVector{<:Real})
    cv = ControlVector(bg_interface.analysis_config)
    @assert length(ξ) == length(cv.values) "Control vector size mismatch with background operator"

    cv.values .= ξ
    result_cv = BackgroundError.apply_background_error_sqrt(bg_interface.background_operator, cv)

    return copy(result_cv.values)
end

"""
    create_gsi_background_error(grid_config, bg_error_config) -> GSIBackgroundErrorInterface

Create realistic GSI background error covariance for atmospheric variables.
"""
function create_gsi_background_error(grid_config, bg_error_config)

    println("Creating GSI background error operator (incremental form)...")

    # Define grid dimensions
    nx, ny, nz = grid_config.nx, grid_config.ny, grid_config.nsig
    n_total = 4 * nx * ny * nz + nx * ny  # u, v, t, q (3D) + ps (2D)

    # Derive analysis configuration used by ControlVector / B-operator
    analysis_config = AnalysisConfig(
        grid_size = (nx, ny, nz),
        ensemble_size = _get_config_value(bg_error_config, :ensemble_size, 0),
        hybrid_coeff = _get_config_value(bg_error_config, :hybrid_coeff, 0.0),
        use_hybrid = _get_config_value(bg_error_config, :use_hybrid, false),
        max_iterations = _get_config_value(bg_error_config, :max_iterations, 200),
        convergence_tol = _get_config_value(bg_error_config, :convergence_tol, 1e-5)
    )

    # Create background error operator using operational factorized form
    berror_file = _get_config_value(bg_error_config, :berror_file, nothing)
    use_balance = _get_config_value(bg_error_config, :use_balance, false)
    use_vertical_modes = _get_config_value(bg_error_config, :use_vertical_modes, true)

    background_operator = BackgroundError.create_operational_bmatrix(
        analysis_config;
        berror_file = berror_file,
        use_balance = use_balance,
        use_vertical_modes = use_vertical_modes
    )

    # Create background error interface
    bg_interface = GSIBackgroundErrorInterface(
        bg_error_config,
        grid_config,
        analysis_config,
        background_operator
    )

    # Create realistic error variances for diagnostics / documentation
    println("  Configuring variance templates (diagnostic only)...")

    t_var = create_temperature_error_variance(nx, ny, nz, bg_error_config)
    u_var = create_wind_error_variance(nx, ny, nz, bg_error_config, :u)
    v_var = create_wind_error_variance(nx, ny, nz, bg_error_config, :v)
    q_var = create_humidity_error_variance(nx, ny, nz, bg_error_config)
    ps_var = create_pressure_error_variance(nx, ny, bg_error_config)

    bg_interface.error_variances[:temperature] = t_var
    bg_interface.error_variances[:u_wind] = u_var
    bg_interface.error_variances[:v_wind] = v_var
    bg_interface.error_variances[:humidity] = q_var
    bg_interface.error_variances[:pressure] = ps_var

    println("  Background error operator configured for state size $(n_total)")
    println("  Balance operators enabled: $(use_balance)")
    println("  Vertical modes enabled: $(use_vertical_modes)")

    return bg_interface
end

"""
    create_temperature_error_variance(nx, ny, nz, config) -> Matrix{Float64}

Create realistic temperature background error variances.
"""
function create_temperature_error_variance(nx, ny, nz, config)
    t_var = zeros(nx, ny, nz)
    
    # Base temperature error (K²)
    base_error = get(config.variance_scaling, :temperature, 2.0)
    
    for k in 1:nz
        # Vertical profile: larger errors in upper atmosphere
        sigma_level = (k - 1) / (nz - 1)
        vertical_factor = 1.0 + 2.0 * sigma_level  # Increase with height
        
        for j in 1:ny, i in 1:nx
            # Horizontal variation (simplified)
            lat_factor = 1.0 + 0.5 * cos(2π * (j - 1) / ny)  # Latitude dependence
            
            t_var[i, j, k] = base_error * vertical_factor * lat_factor
        end
    end
    
    return t_var
end

"""
    create_wind_error_variance(nx, ny, nz, config, component) -> Matrix{Float64}

Create realistic wind component background error variances.
"""
function create_wind_error_variance(nx, ny, nz, config, component)
    wind_var = zeros(nx, ny, nz)
    
    # Base wind error (m²/s²)
    base_error = get(config.variance_scaling, :wind, 1.5)
    
    for k in 1:nz
        # Vertical profile: jet stream levels have higher errors
        sigma_level = (k - 1) / (nz - 1)
        if sigma_level > 0.6 && sigma_level < 0.9  # Jet stream region
            vertical_factor = 2.0
        else
            vertical_factor = 1.0
        end
        
        for j in 1:ny, i in 1:nx
            wind_var[i, j, k] = base_error * vertical_factor
        end
    end
    
    return wind_var
end

"""
    create_humidity_error_variance(nx, ny, nz, config) -> Matrix{Float64}

Create realistic humidity background error variances.
"""
function create_humidity_error_variance(nx, ny, nz, config)
    q_var = zeros(nx, ny, nz)
    
    # Base humidity error (kg²/kg²)
    base_error = get(config.variance_scaling, :humidity, 0.8)
    
    for k in 1:nz
        # Vertical profile: higher errors in boundary layer
        sigma_level = (k - 1) / (nz - 1)
        if sigma_level < 0.3  # Boundary layer
            vertical_factor = 3.0
        else
            vertical_factor = 1.0
        end
        
        for j in 1:ny, i in 1:nx
            q_var[i, j, k] = base_error * vertical_factor
        end
    end
    
    return q_var
end

"""
    create_pressure_error_variance(nx, ny, config) -> Matrix{Float64}

Create realistic surface pressure background error variances.
"""
function create_pressure_error_variance(nx, ny, config)
    ps_var = zeros(nx, ny)
    
    # Base pressure error (Pa²)
    base_error = get(config.variance_scaling, :pressure, 1.2)
    
    for j in 1:ny, i in 1:nx
        # Higher errors over land/complex terrain (simplified)
        terrain_factor = 1.0 + 0.5 * rand()  # Randomized terrain effect
        ps_var[i, j] = base_error * terrain_factor
    end
    
    return ps_var
end

"""
    create_correlation_matrix(bg_interface, nx, ny, nz) -> SparseMatrixCSC

Create sparse correlation matrix with realistic atmospheric correlation structures.
"""
function create_correlation_matrix(bg_interface, nx, ny, nz)
    
    n_vars = 5  # u, v, t, q, ps
    n_total = nx * ny * nz * (n_vars - 1) + nx * ny  # 3D vars + 2D pressure
    
    # Use sparse matrix construction for efficiency
    I_indices = Int[]
    J_indices = Int[]
    values = Float64[]
    
    println("    Creating correlation structure for $(n_total) state variables...")
    
    # Horizontal correlation length (grid points)
    h_corr = bg_interface.correlation_lengths[:horizontal] / 12.0  # 12 km grid spacing
    v_corr = bg_interface.correlation_lengths[:vertical]
    
    # Build correlations (simplified Gaussian model)
    for var in 1:(n_vars-1)  # 3D variables
        var_offset = (var - 1) * nx * ny * nz
        
        for k in 1:nz, j in 1:ny, i in 1:nx
            idx = var_offset + (k-1)*nx*ny + (j-1)*nx + i
            
            # Self correlation
            push!(I_indices, idx)
            push!(J_indices, idx)
            push!(values, 1.0)
            
            # Horizontal correlations
            for dj in -2:2, di in -2:2
                if di == 0 && dj == 0
                    continue
                end
                
                i2, j2 = i + di, j + dj
                if 1 <= i2 <= nx && 1 <= j2 <= ny
                    idx2 = var_offset + (k-1)*nx*ny + (j2-1)*nx + i2
                    
                    dist = sqrt(di^2 + dj^2)
                    correlation = exp(-0.5 * (dist / h_corr)^2)
                    
                    if correlation > 0.1  # Threshold for sparsity
                        push!(I_indices, idx)
                        push!(J_indices, idx2)
                        push!(values, correlation)
                    end
                end
            end
            
            # Vertical correlations
            for dk in [-1, 1]
                k2 = k + dk
                if 1 <= k2 <= nz
                    idx2 = var_offset + (k2-1)*nx*ny + (j-1)*nx + i
                    
                    correlation = exp(-0.5 * (abs(dk) / v_corr)^2)
                    
                    if correlation > 0.1
                        push!(I_indices, idx)
                        push!(J_indices, idx2)
                        push!(values, correlation)
                    end
                end
            end
        end
    end
    
    # Surface pressure correlations (2D field)
    ps_offset = (n_vars - 1) * nx * ny * nz
    for j in 1:ny, i in 1:nx
        idx = ps_offset + (j-1)*nx + i
        
        # Self correlation
        push!(I_indices, idx)
        push!(J_indices, idx)
        push!(values, 1.0)
        
        # Horizontal correlations
        for dj in -2:2, di in -2:2
            if di == 0 && dj == 0
                continue
            end
            
            i2, j2 = i + di, j + dj
            if 1 <= i2 <= nx && 1 <= j2 <= ny
                idx2 = ps_offset + (j2-1)*nx + i2
                
                dist = sqrt(di^2 + dj^2)
                correlation = exp(-0.5 * (dist / h_corr)^2)
                
                if correlation > 0.1
                    push!(I_indices, idx)
                    push!(J_indices, idx2)
                    push!(values, correlation)
                end
            end
        end
    end
    
    # Create sparse matrix
    B_matrix = sparse(I_indices, J_indices, values, n_total, n_total)
    
    println("    Correlation matrix created with $(nnz(B_matrix)) non-zero elements")
    
    return B_matrix
end

# =============================================================================
# GSI Model Interface
# =============================================================================

"""
    GSIModelInterface

Interface for integrating atmospheric forecast models with DRP-4DVar.
"""
struct GSIModelInterface
    model_type::Symbol  # :wrf, :gfs, :fv3, etc.
    time_step::Float64  # Model time step (hours)
    
    # Model operators by time step
    tangent_linear_ops::Dict{Int, Any}
    adjoint_ops::Dict{Int, Any}
    nonlinear_ops::Dict{Int, Any}
    
    function GSIModelInterface(model_type::Symbol, time_step::Float64=1.0)
        new(model_type, time_step, 
            Dict{Int, Any}(), Dict{Int, Any}(), Dict{Int, Any}())
    end
end

"""
    add_model_operators!(interface, time_steps, operator_type)

Add model operators for specified time steps.
"""
function add_model_operators!(interface::GSIModelInterface, 
                             time_steps::Vector{Int},
                             operator_type::Symbol=:linear)
    
    for time_idx in time_steps
        if operator_type == :linear || operator_type == :all
            # Create simplified tangent linear operator
            interface.tangent_linear_ops[time_idx] = create_linear_model_operator(
                interface.model_type, time_idx, interface.time_step
            )
        end
        
        if operator_type == :nonlinear || operator_type == :all
            # Create nonlinear model operator
            interface.nonlinear_ops[time_idx] = create_nonlinear_model_operator(
                interface.model_type, time_idx, interface.time_step
            )
        end
    end
    
    return nothing
end

"""
    create_linear_model_operator(model_type, time_idx, dt)

Create linearized atmospheric forecast model operator.
"""
function create_linear_model_operator(model_type::Symbol, time_idx::Int, dt::Float64)
    
    # Simplified linear model operator (identity + small perturbation)
    # In practice, this would be the tangent linear of WRF/GFS
    
    return function(state_vector::Vector{Float64})
        n = length(state_vector)
        
        # Simple advection-diffusion-like operator
        A = Matrix{Float64}(I, n, n)
        
        # Add some dynamics (very simplified)
        damping = 0.95  # Slight damping
        for i in 1:n
            A[i, i] = damping
        end
        
        # Add coupling between variables (simplified)
        n_per_field = n ÷ 5
        for i in 1:min(n_per_field, n-n_per_field)
            # Temperature-wind coupling
            A[i, 2*n_per_field + i] = 0.01  # T affects u
            A[n_per_field + i, 2*n_per_field + i] = 0.01  # T affects v
        end
        
        return A * state_vector
    end
end

"""
    create_nonlinear_model_operator(model_type, time_idx, dt)

Create nonlinear atmospheric forecast model operator.
"""
function create_nonlinear_model_operator(model_type::Symbol, time_idx::Int, dt::Float64)
    
    return function(state_vector::Vector{Float64})
        # Apply linear operator first
        linear_op = create_linear_model_operator(model_type, time_idx, dt)
        linear_result = linear_op(state_vector)
        
        # Add nonlinear terms (very simplified)
        nonlinear_correction = 0.01 * sin.(0.1 * linear_result)
        
        return linear_result + nonlinear_correction
    end
end

# =============================================================================
# Atmospheric Test Case Generator
# =============================================================================

"""
    AtmosphericTestCase

Container for realistic atmospheric data assimilation test cases.
"""
struct AtmosphericTestCase
    name::String
    grid_config::Any
    background_field::Vector{Float64}
    background_error::GSIBackgroundErrorInterface
    observations::Dict{Int, Vector{Float64}}
    observation_operators::Dict{Int, Any}
    model_operators::Dict{Int, Any}
    observation_errors::Dict{Int, Matrix{Float64}}
    
    # Metadata
    domain_type::Symbol  # :regional, :global
    resolution::String
    variables::Vector{Symbol}
    time_window::Int
    n_observations::Int
end

"""
    create_atmospheric_test_case(case_name; options...)

Create realistic atmospheric data assimilation test case.
"""
function create_atmospheric_test_case(case_name::String="regional_conus";
                                    nx::Int=100, ny::Int=80, nz::Int=30,
                                    time_window::Int=6,
                                    n_obs_per_time::Int=500,
                                    obs_types::Vector{Symbol}=[:temperature, :wind, :pressure])
    
    println("Creating atmospheric test case: $case_name")
    
    # Create grid configuration
    grid_config = create_test_grid_config(nx, ny, nz, case_name)
    
    # Create realistic background field
    background_field = create_test_background_field(nx, ny, nz)
    
    # Create realistic background error covariance
    bg_error_config = BackgroundError.BackgroundErrorConfig(
        variance_scaling = Dict(
            :temperature => 4.0,  # K²
            :wind => 2.25,        # (m/s)²
            :humidity => 1.0,     # (kg/kg)²  
            :pressure => 150.0    # Pa²
        ),
        correlation_lengths = Dict(
            :horizontal => 120.0,  # km
            :vertical => 2.0       # scale heights
        )
    )
    
    bg_error_interface = create_gsi_background_error(grid_config, bg_error_config)
    
    # Create observation interface
    obs_interface = GSIObservationInterface(grid_config)
    
    # Add observations for each time step
    observations = Dict{Int, Vector{Float64}}()
    obs_errors = Dict{Int, Matrix{Float64}}()
    
    for time_idx in 1:time_window
        
        # Generate observation locations
        obs_locations = generate_observation_locations(nx, ny, nz, n_obs_per_time, case_name)
        
        # Add different observation types
        if :temperature in obs_types
            add_temperature_observations!(obs_interface, obs_locations, [time_idx])
        end
        
        if :wind in obs_types
            add_wind_observations!(obs_interface, obs_locations, [time_idx], :u)
            add_wind_observations!(obs_interface, obs_locations, [time_idx], :v)
        end

        if :pressure in obs_types
            add_pressure_observations!(obs_interface, obs_locations, [time_idx])
        end

        # Create synthetic observations
        observations[time_idx] = create_synthetic_observations(
            background_field, obs_locations, obs_types, time_idx, grid_config
        )
        
        # Create observation error covariance
        n_obs = length(observations[time_idx])
        obs_errors[time_idx] = create_observation_error_matrix(n_obs, obs_types)
    end
    
    # Create observation and model operators
    observation_operators = create_observation_operators(obs_interface)
    model_interface = GSIModelInterface(:wrf, 1.0)  # 1-hour time step
    add_model_operators!(model_interface, collect(1:time_window), :all)
    model_operators = model_interface.tangent_linear_ops
    
    # Create test case
    test_case = AtmosphericTestCase(
        case_name,
        grid_config,
        background_field,
        bg_error_interface,
        observations,
        observation_operators,
        model_operators,
        obs_errors,
        :regional,
        "$(grid_config.dx/1000)km",
        [:u, :v, :t, :q, :ps],
        time_window,
        sum(length(obs) for obs in values(observations))
    )
    
    println("Test case created successfully:")
    println("  Domain: $(test_case.domain_type) $(test_case.resolution)")
    println("  Grid: $(nx)×$(ny)×$(nz)")
    println("  Time window: $(time_window) hours")
    println("  Total observations: $(test_case.n_observations)")
    
    return test_case
end

"""
    create_test_grid_config(nx, ny, nz, case_name)

Create grid configuration for test cases.
"""
function create_test_grid_config(nx, ny, nz, case_name)
    if case_name == "regional_conus"
        # CONUS domain configuration
        return (
            nx = nx, ny = ny, nsig = nz,
            dx = 12000.0, dy = 12000.0,  # 12 km resolution
            center_lon = -95.0, center_lat = 40.0,
            proj = "lambert_conformal"
        )
    elseif case_name == "regional_europe"
        return (
            nx = nx, ny = ny, nsig = nz,
            dx = 10000.0, dy = 10000.0,  # 10 km resolution  
            center_lon = 10.0, center_lat = 50.0,
            proj = "lambert_conformal"
        )
    else
        # Default regional domain
        return (
            nx = nx, ny = ny, nsig = nz,
            dx = 15000.0, dy = 15000.0,  # 15 km resolution
            center_lon = 0.0, center_lat = 45.0,
            proj = "lambert_conformal"
        )
    end
end

"""
    create_test_background_field(nx, ny, nz) -> Vector{Float64}

Create realistic background atmospheric field.
"""
function create_test_background_field(nx, ny, nz)
    
    # Total state size: u, v, t, q (3D) + ps (2D)
    n_3d = nx * ny * nz
    n_2d = nx * ny
    n_total = 4 * n_3d + n_2d
    
    background = zeros(n_total)
    
    println("  Creating background field with $(n_total) state variables...")
    
    # U-wind component (m/s)
    u_field = reshape(background[1:n_3d], nx, ny, nz)
    for k in 1:nz, j in 1:ny, i in 1:nx
        # Jet stream profile + random perturbations
        sigma = (k - 1) / (nz - 1)
        lat_factor = sin(π * (j - 1) / (ny - 1))
        u_field[i, j, k] = 20.0 * lat_factor * exp(-((sigma - 0.7)/0.2)^2) + randn() * 2.0
    end
    background[1:n_3d] = vec(u_field)
    
    # V-wind component (m/s) - smaller magnitude
    v_field = reshape(background[(n_3d+1):(2*n_3d)], nx, ny, nz)
    for k in 1:nz, j in 1:ny, i in 1:nx
        v_field[i, j, k] = randn() * 5.0
    end
    background[(n_3d+1):(2*n_3d)] = vec(v_field)
    
    # Temperature (K)
    t_field = reshape(background[(2*n_3d+1):(3*n_3d)], nx, ny, nz)
    for k in 1:nz, j in 1:ny, i in 1:nx
        # Standard atmosphere + perturbations
        sigma = (k - 1) / (nz - 1)
        t_field[i, j, k] = 288.15 - 65.0 * sigma + randn() * 3.0  # Lapse rate
    end
    background[(2*n_3d+1):(3*n_3d)] = vec(t_field)
    
    # Specific humidity (kg/kg)
    q_field = reshape(background[(3*n_3d+1):(4*n_3d)], nx, ny, nz)
    for k in 1:nz, j in 1:ny, i in 1:nx
        # Exponential decrease with height
        sigma = (k - 1) / (nz - 1)
        q_field[i, j, k] = 0.015 * exp(-3.0 * sigma) + abs(randn()) * 0.002
    end
    background[(3*n_3d+1):(4*n_3d)] = vec(q_field)
    
    # Surface pressure (Pa)
    ps_field = reshape(background[(4*n_3d+1):end], nx, ny)
    for j in 1:ny, i in 1:nx
        ps_field[i, j] = 101325.0 + randn() * 500.0  # Sea level pressure + noise
    end
    background[(4*n_3d+1):end] = vec(ps_field)
    
    println("  Background field statistics:")
    println("    U-wind: $(mean(background[1:n_3d])) ± $(std(background[1:n_3d])) m/s")
    println("    V-wind: $(mean(background[(n_3d+1):(2*n_3d)])) ± $(std(background[(n_3d+1):(2*n_3d)])) m/s")
    println("    Temperature: $(mean(background[(2*n_3d+1):(3*n_3d)])) ± $(std(background[(2*n_3d+1):(3*n_3d)])) K")
    println("    Humidity: $(mean(background[(3*n_3d+1):(4*n_3d)])) ± $(std(background[(3*n_3d+1):(4*n_3d)])) kg/kg")
    println("    Surface pressure: $(mean(background[(4*n_3d+1):end])) ± $(std(background[(4*n_3d+1):end])) Pa")
    
    return background
end

"""
    generate_observation_locations(nx, ny, nz, n_obs, case_name) -> Matrix{Float64}

Generate realistic observation locations for the test case.
"""
function generate_observation_locations(nx, ny, nz, n_obs, case_name)
    
    locations = zeros(n_obs, 3)  # [i, j, k] grid coordinates
    
    for obs_idx in 1:n_obs
        # Random horizontal location
        locations[obs_idx, 1] = 1 + rand() * (nx - 1)  # i coordinate
        locations[obs_idx, 2] = 1 + rand() * (ny - 1)  # j coordinate
        
        # Vertical coordinate (pressure level)
        # Most observations in lower atmosphere
        if rand() < 0.6  # 60% surface/boundary layer
            locations[obs_idx, 3] = 1 + rand() * (nz * 0.3)
        elseif rand() < 0.8  # 20% mid-troposphere
            locations[obs_idx, 3] = nz * 0.3 + rand() * (nz * 0.4)
        else  # 20% upper troposphere
            locations[obs_idx, 3] = nz * 0.7 + rand() * (nz * 0.3)
        end
    end
    
    return locations
end

"""
    create_synthetic_observations(background, locations, obs_types, time_idx) -> Vector{Float64}

Create synthetic observations by applying operators to background + truth perturbations.
"""
function create_synthetic_observations(background, locations, obs_types, time_idx, grid)
    
    n_obs_total = size(locations, 1) * length(obs_types)
    observations = zeros(n_obs_total)
    
    obs_idx = 1
    
    for obs_type in obs_types
        for loc_idx in 1:size(locations, 1)
            i, j, k = round.(Int, locations[loc_idx, :])
            nx, ny, nz = _grid_dimensions(grid)
            i = clamp(i, 1, nx)
            j = clamp(j, 1, ny)
            k = clamp(k, 1, nz)
            
            if obs_type == :temperature
                # Extract temperature value with interpolation (simplified)
                observations[obs_idx] = extract_field_value(background, :t, i, j, k, grid) + randn() * 2.0
            elseif obs_type == :wind
                # Wind speed magnitude
                u_val = extract_field_value(background, :u, i, j, k, grid)
                v_val = extract_field_value(background, :v, i, j, k, grid)
                observations[obs_idx] = sqrt(u_val^2 + v_val^2) + abs(randn()) * 1.5
            elseif obs_type == :pressure
                # Surface pressure
                observations[obs_idx] = extract_field_value(background, :ps, i, j, 1, grid) + randn() * 100.0
            end
            
            obs_idx += 1
        end
    end
    
    return observations
end

"""
    extract_field_value(background, field_type, i, j, k) -> Float64

Extract field value from background state vector.
"""
function extract_field_value(background, field_type, i, j, k, grid)
    nx, ny, nz = _grid_dimensions(grid)
    n_3d = nx * ny * nz
    
    i = clamp(i, 1, nx)
    j = clamp(j, 1, ny)
    k = clamp(k, 1, nz)

    if field_type == :u
        field_offset = 0
    elseif field_type == :v
        field_offset = n_3d
    elseif field_type == :t
        field_offset = 2 * n_3d
    elseif field_type == :q
        field_offset = 3 * n_3d
    elseif field_type == :ps
        # Surface pressure (2D field)
        ps_offset = 4 * n_3d
        ps_idx = ps_offset + (j-1)*nx + i
        return ps_idx <= length(background) ? background[ps_idx] : 0.0
    else
        return 0.0
    end
    
    # 3D field index
    field_idx = field_offset + (k-1)*nx*ny + (j-1)*nx + i
    
    # Bounds checking
    if 1 <= field_idx <= length(background)
        return background[field_idx]
    else
        return 0.0
    end
end

"""
    create_observation_error_matrix(n_obs, obs_types) -> Matrix{Float64}

Create diagonal observation error covariance matrix.
"""
function create_observation_error_matrix(n_obs, obs_types)
    
    R = zeros(n_obs, n_obs)
    
    obs_idx = 1
    n_locs = n_obs ÷ length(obs_types)
    
    for obs_type in obs_types
        for loc_idx in 1:n_locs
            if obs_type == :temperature
                R[obs_idx, obs_idx] = 4.0  # 2 K standard error
            elseif obs_type == :wind
                R[obs_idx, obs_idx] = 2.25  # 1.5 m/s standard error
            elseif obs_type == :pressure
                R[obs_idx, obs_idx] = 10000.0  # 100 Pa standard error
            else
                R[obs_idx, obs_idx] = 1.0  # Default
            end
            obs_idx += 1
        end
    end
    
    return R
end

# =============================================================================
# Operational DRP-4DVar Interface
# =============================================================================

"""
    run_operational_drp4dvar(test_case; ensemble_size=40, optimizer="lbfgs")

Run DRP-4DVar analysis on a realistic atmospheric test case.
"""
function run_operational_drp4dvar(test_case::AtmosphericTestCase;
                                 ensemble_size::Int=40,
                                 optimizer::String="lbfgs",
                                 max_outer_loops::Int=3,
                                 convergence_tolerance::Float64=1e-6,
                                 use_localization::Bool=true,
                                 localization_radius::Float64=200.0,
                                 max_inner_loops::Int=150,
                                 ensemble_inflation::Float64=1.02,
                                 adaptive_inflation::Bool=true)
    
    println("="^80)
    println("OPERATIONAL DRP-4DVAR ANALYSIS")
    println("Test Case: $(test_case.name)")
    println("="^80)
    
    # Initialize DRP-4DVar with operational settings
    drp4dvar = DRP4DVar(
        ensemble_size = ensemble_size,
        max_outer_loops = max_outer_loops,
        max_inner_loops = max_inner_loops,
        convergence_tolerance = convergence_tolerance,
        time_window = test_case.time_window,
        optimizer = optimizer,
        use_localization = use_localization,
        localization_radius = localization_radius,
        ensemble_inflation = ensemble_inflation,
        adaptive_inflation = adaptive_inflation
    )
    
    # Run DRP-4DVar analysis
    analysis_state, statistics = run_drp4dvar(
        drp4dvar,
        test_case.background_field,
        test_case.background_error,
        test_case.observations,
        test_case.observation_operators,
        test_case.model_operators
    )
    
    # Compute analysis diagnostics
    analysis_increment = analysis_state - test_case.background_field
    increment_norm = norm(analysis_increment)
    rms_increment = sqrt(mean(analysis_increment.^2))
    
    # Innovation statistics
    innovation_stats = compute_innovation_statistics(
        test_case, analysis_state, test_case.observation_operators
    )
    
    println("\nANALYSIS RESULTS:")
    println("  Final cost function: $(get(statistics, "final_cost", NaN))")
    println("  Analysis increment norm: $(increment_norm)")
    println("  RMS analysis increment: $(rms_increment)")
    println("  Execution time: $(get(statistics, "total_execution_time", NaN)) seconds")
    println("  Convergence achieved: $(get(statistics, "converged", false))")
    
    println("\nINNOVATION STATISTICS:")
    for (key, value) in innovation_stats
        println("  $(key): $(value)")
    end
    
    # Return comprehensive results
    operational_results = Dict(
        "analysis_state" => analysis_state,
        "analysis_increment" => analysis_increment,
        "background_state" => test_case.background_field,
        "statistics" => statistics,
        "innovation_statistics" => innovation_stats,
        "test_case_info" => Dict(
            "name" => test_case.name,
            "domain_type" => test_case.domain_type,
            "resolution" => test_case.resolution,
            "n_observations" => test_case.n_observations,
            "time_window" => test_case.time_window
        ),
        "method_parameters" => Dict(
            "ensemble_size" => ensemble_size,
            "optimizer" => optimizer,
            "use_localization" => use_localization,
            "localization_radius" => localization_radius
        )
    )
    
    println("\nOperational DRP-4DVar analysis completed successfully!")
    println("="^80)
    
    return operational_results
end

"""
    compute_innovation_statistics(test_case, analysis_state, obs_operators)

Compute observation-minus-forecast statistics for analysis validation.
"""
function compute_innovation_statistics(test_case, analysis_state, obs_operators)
    
    innovation_stats = Dict{String, Float64}()
    
    total_obs = 0
    total_innovation_sq = 0.0
    
    for (time_idx, observations) in test_case.observations
        if haskey(obs_operators, time_idx)
            # Compute model equivalent for analysis
            h_analysis = obs_operators[time_idx](analysis_state)
            
            # Compute innovations (obs - analysis)
            innovations = observations - h_analysis
            
            total_obs += length(innovations)
            total_innovation_sq += sum(innovations.^2)
        end
    end
    
    if total_obs > 0
        innovation_stats["total_observations"] = Float64(total_obs)
        innovation_stats["rms_innovation"] = sqrt(total_innovation_sq / total_obs)
        innovation_stats["mean_abs_innovation"] = sqrt(total_innovation_sq / total_obs)  # Simplified
    end
    
    return innovation_stats
end

end # module GSIIntegration
